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Abstract. We show in this paper the results on the phase transition of the so-called 
fully frustrated simple cubic lattice with the Ising spin model. We use here the Monte 
Carlo method with the flat energy-histogram Wang-Landau technique which is very 
powerful to detect weak first-order phase transition. We show that the phase transition 
is clearly of first order, providing a definite answer to a question raised 25 years ago. 
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1. Introduction 

Statistical physics provides powerful methods to study behaviors of systems of 
interacting particles. In particular, different kinds of transition from one phase to 
another has been studied with efficiency during the last 40 years by exact methods pQ, 
renormalization group, high- low-temperature expansions [2], numerical simulations, ... 
Experiments have verified most of these theoretical results. Among the most studied 
subjects, we mention the effect of the frustration in spin systems. The frustration 
is known to be the origin of spectacular properties such as large ground state (GS) 
degeneracy, successive phase transitions, partially disordered phase, reentrance and 
disorder lines. Though these aspects have been found in exactly solved models[3], we 
believe that many of these features remain in complicated frustrated systems where 
exact solutions are not available. These general frustrated systems still constitute at 
present a challenge for theoretical physics [I]. 

Let us recall the definition of a frustrated system. When a spin cannot fully satisfy 
energetically all the interactions with its neighbors, it is "frustrated". This occurs when 
the interactions are in competition with each other, for instance incompatible nearest- 
neighbor (NN) and next-nearest-neighbor (NNN) interactions, or when the lattice 
geometry does not allow a spin to satisfy all interaction bonds simultaneously such as 
the triangular antiferromagnet. Except a few two-dimensional frustrated Ising systems 
where exact methods have been devised to solve with mathematical elegance [H El EJ [7J 
El E] , most systems have recourse to numerical simulations and various approximations. 
One of the most studied systems is the stacked triangular antiferromagnet (STA) with 
interaction between NN. This system with Ising[10j, XY and Heisenberg spins [12] 
have been intensively studied since 1987(131 El [151 OS HO US CEHl ED] , but only recently 
that the 20-year controversy comes to an end [211 [22l [231 121 [251 [261 1211 M\ [291 [30] . Note 
that numerical simulations require now new efficient algorithms to deal with frustrated 
systems [23 [30]. 

There is another fully frustrated system. Initially defined in two-dimensions (2D) 
on a square lattice by Villain [31], this model has been generalized in three dimensions 
(3D) as shown in Fig. [I] by Blankschtein et al. [32]. A detailed description of the model 
will be presented in section EJ The nature of the phase transition in the classical XY[33J 
and Heisenberg [34J spin models has been recently investigated. It was shown that it is 
a first-order transition putting an end to a 25- year long controversial issue |35[ 136] . In 
this paper, we extend our study to the case of Ising spin model. 

In Section [2] we describe the model and give some technical details of the Wang- 
Landau (WL) methods as applied in the present paper. Section [3] shows our results. 
Concluding remarks are given in section [H 
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2. Model and Wang-Landau Method 

The model shown in Fig. [1] has been previously called "fully frustrated simple cubic 
lattice" (FFSCL) by Blankschtein et al.|32]. The Hamiltonian is given by 

W = ~ / ] JijSi-Sj, (1) 

(vi) 

where Sj is the Ising spin of values ±1 at the lattice site i, £V ^ is made over the NN spin 
pairs Si and Sj with interaction J^. We take J^- = — J (J > 0) for antiferromagnetic 
bonds indicated by discontinued lines in Fig. [H and Jy = J for ferromagnetic bonds 
indicated by continued lines. The 2D Villain's model has been intensively studied 
with Ising model[3U |37] and XY spin model due to its application in arrays of planar 
Josephson's j unctions [39| HUl HI]. 

Let us recall some results on the present model. The GS degeneracy is infinite due 
to the fact that each face of the cube is frustrated, there is thus an infinite number to 
arrange the spins in an infinite crystal. Note that ferromagnetic state is one of the GS 
spin configurations. In an early MC study [12], it has been shown that as the temperature 
T increases, the system selects the long-range ferromagnetic state at low T but goes to 
a partially disordered phase where two of the 8 sublattices of the cube are disordered. 
The passage to this phase does not have the characteristics of a phase transition. The 
specific heat shows a "shoulder" at T ~ 0.5 (in unit of J/ks), far below the transition 
temperature for the whole system occurring at T ~ 1.345. Note that the nature of 
the low-T ordering of the present model is still not elucidated. However, in 1987 we 
have shown [5] in an exactly solved 2D model that a partial disorder can coexist with an 
order at equilibrium. Therefore, we believe that the early observation of two disordered 
sublattices in an ordered phase may have the same origin rooted in the frustration and 
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in the order selection by entropy j2D EH EH]- The shoulder of the specific heat may 
turn out to be a true phase transition. This point has to be checked with careful MC 
simulations using very large lattice sizes. This is a formidable task, but it is not the 
purpose of this work. In the present work, we concentrate our attention on the nature 
of the overall phase transition occurring at a higher temperature. Using the Landau- 
Ginzburg- Wilson theory, Blankschtein et al.[32]have found a weak first-order transition. 
Our previous work in 1985 using a standard MC algorithm with short runs and small 
lattice sizes permitted by the computer capacity at that time [42] show a second-order 
transition with an unusual critical properties in contradiction with the prediction of 
Blankschtein et al. In the light of new results on frustrated systems obtained not only by 
new efficient MC algorithms but also by today's huge computer capacity [29] |30 | |33 | 134"] . 
we study this problem again in order to get a definite answer to that question. 

For weak first-order transitions, MC simulations with the standard Metropolis 
algorithm cannot give results with good precision even with the use of large sizes and 
long runs. This is because the algorithm does not allow us, among other difficulties, to 
easily sample rare microscopic states. Wang and Landau [43] have recently proposed a 
MC algorithm which allowed to study classical statistical models with difficultly accessed 
microscopic states. In particular, it permits to detect with efficiency weak first-order 
transitions [29 [ 130] [33] The algorithm uses a random walk in energy space in order to 
obtained an accurate estimate for the density of states g(E) which is defined as the 
number of spin configurations for any given E. This method is based on the fact that 
a flat energy histogram H(E) is produced if the probability for the transition to a state 
of energy E is proportional to g(E)^ 1 . At the beginning of the simulation, the density 
of states (DOS) is set equal to one for all possible energies, g(E) = 1. We begin a 
random walk in energy space (E) by choosing a site randomly and flipping its spin with 
a probability proportional to the inverse of the temporary density of states (DOS). In 
general, if E and E' are the energies before and after a spin is flipped, the transition 
probability from E to E' is 



Of course, to enhance the possibility to access to rare states, some tricks have been 
devised. Each time an energy level E is visited, the DOS is modified by a modification 
factor / > whether the spin flipped or not, i.e. g(E) — >■ g(E)f. At the beginning of the 
random walk, the modification factor / can be as large as e 1 ~ 2.7182818. A histogram 
H{E) records the number of times a state of energy E is visited. Each time the energy 
histogram satisfies a certain "flatness" criterion, / is reduced according to / — > \ff and 
H(E) is reset to zero for all energies. The reduction process of the modification factor 
/ is repeated several times until a final value /fi na i which close enough to one. The 
histogram is considered as flat if 



for all energies, where x% is chosen between 70% and 95% and (H(E)) is the average 
histogram. 



p(E^E') = mm[g(E)/g(E'),l}. 



(2) 



H(E)>x%.(H(E)) 



(3) 
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The WL method has been applied to our spin models with success in our recent 
papers |29j [301 We emphasize that for efficiency, we consider here a multi subinterval 
energy scale within an energy range of interest [4"5l l4"6] (E min , £ max ) which covers not all 
possible energies of the system but all energies in the region will will use in applications. 
We divide this energy range to R subintervals, the minimum energy of the i — th 
subinterval is -E^in (* = 1,2,..., i?), and the maximum is E l m&x = E^t, + 2AE, where 
AE can be chosen large enough for a smooth boundary between two subintervals. The 
WL algorithm is used to calculate the relative DOS of each subinterval (B^, -E^ax) 
with a flatness criterion x% = 95%. Note that we reject a spin flip and do not update 
g(E) and the energy histogram H(E) of the current energy level E if the spin-flip trial 
would result in an energy outside the energy segment. The DOS of the whole range is 
obtained by joining the DOS of each subinterval (E^ n + AE, -E^ax — AE). 

The thermodynamic quant it ies|43[ Hlj can be evaluated by 

(E n ) = \Y. En 9( E ) eM~E/k B T) (4) 

C ~ ^~^ 2 (5) 
Cv ~ k B T 2 (5) 

(M n ) = \Y. Mn 9( E ) zM-E/k B T) (6) 
(M 2 ) - (M) 2 

X = ksT (7) 

where Z is the partition function defined by Z = g{E) exp(—E/kBT). The 
canonical distribution at a temperature T can be calculated simply by P(E, T) = 
±g{E)exp{-E/k B T). 

The simulations have been carried our on a rack of several hundreds of 64-bit CPU. 
For a given size L, the calculation takes, depending on L, from a few weeks to several 
months to have the required histogram flatness. 



3. Results 



We have started the simulations from the system linear size L = 60 (the system size is 
I/ 3 ). But only from L = 90 that a sign of first-order transition appears. Therefore, we 
use extremely large sizes up to 180. Periodic boundary conditions are used in the three 
directions. J = 1 is taken as the unit of energy in the following. 

We show in Fig. [2] the energy per spin and the specific heat, for L = 180, using the 
flat histogram obtained with WL method. Several remarks are in order: 

i) the energy at the largest size shows a 'pseudo" discontinuity at the transition 
temperature Tc — 1.34814. We will see below that this discontinuity is confirmed by 
the double-peak energy histogram at this temperature, 

ii) the specific heat shows a very strong size dependence. It should be noted that 
the specific heat is calculated from the fluctuation of the energy of the system at a 
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given T [see Eq. (151)]. not by the derivative of E with respect to T. Therefore, when 
the energy has a discontinuity at Tq with two levels E\ and E 2 , the average energy 
is E = (Ei + E-i)j2. It is the fluctuations of E which gives rise to Cy, and we will 
not see a delta-like function should Cy is calculated by the energy derivative. This 
is the reason why in standard MC simulations with the Metropolis algorithm, we do 
not see discontinuity in energy for weak first-order transition (what is sorted out of the 
simulation is an average energy). So, an energy histogram is really needed if we want 
to see weak first order. 

The energy histogram can be realized directly in the old fashion manner by 
measuring the system energy at a given T[47]. However, when relaxation between rare 
states are very slow, we need the temperature-independent WL flat histogram technique 
as described above. We show the WL result in Fig. [3] As seen, for L = 120, the energy 
histogram begins to show a sign of the double-peak structure. The dip between the two 
maxima becomes deeper with increasing size. Note that a "true" discontinuity happens 
only when the dip comes down to E — 0. This requires sizes much larger than L = 180. 
But for our present purpose, we need not to study sizes larger than L = 180. 

We note that the distance between the two peaks, i. e. the latent heat, increases 
with increasing size and reaches ~ 0.005 for L = 180. This is very small compared to 
the value ~ 0.03 for the XY case at L — 48, and to ~ 0.0085 for the Heisenberg case 
at L = 90. The smallness of the latent heat in the present Ising case explains why one 
should go to an extremely large lattice size to detect the first-order transition. 

Let us show in Fig. H] the maximum of Cy versus L in a In — In scale, we find a 
straight line within statistical errors (by a mean least-square fit) with a slope equal to 
d> = 2.794 ± 0.198. This means that C^ ax = AL^ where A is a constant and very 
close to the system dimension d = 3. The fact that Cy ax is proportional to the system 
volume gives another strong signature of a first-order transition. 

The weak first-order transition found here is thus in agreement with the Landau- 
Ginzburg- Wilson theory [32]. This is rather surprising because in other frustrated 
systems such as the STA mentioned in the Introduction, the renormalization group 
with low-order developments in e did not work properly. 

4. Concluding Remarks 

We have showed in this paper the results obtained by the WL flat energy-histogram 
technique on the phase transition in the Ising fully frustrated simple cubic lattice. We 
found that the transition is clearly of first order. Note that the first-order character 
is so weak that it has been observed only at extremely large lattice sizes. This finding 
shows that early studies using standard MC algorithm with short runs and much smaller 
sizes [42J are not correct. Our result confirms the prediction by the Landau- Ginzburg- 
Wilson analysis [32J putting an end to an uncertainty which has lasted for 25 years. 
Together with our recent results [33], [34], we conclude that the fully frustrated simple 
cubic lattice undergoes a first-order transition for Ising, XY and Heisenberg spin models. 
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1.3476 1.348 1.3484 1.3488 T 

Figure 2. Energy per spin E versus temperature T at the lattice size 180 3 (upper 
figure) and specific heat per spin Cy versus T for lattice sizes 120 3 , 140 3 , 160 3 , 180 3 
(lower figure). See text for comments. 
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Figure 4. Maximum of the specific heat C™" 1 versus L in the In — In scale. The 
straight line is a mean least square fit. The slope is <\> = 2.794(i98). Note that the 
specific heat shown in Fig. [5] has been calculated from the fluctuations of the energy. 



It is worth to mention that several other frustrated systems also show a first-order 
transition such as helimagnets^S], FCC[l9] and HCPjHn] antiferromagnets. 

This study shows that one has to be very careful in studying complex systems by 
MC simulations: in some cases such as the one studied here, sizes as large as 80 3 are 
still not sufficient to get a correct conclusion. Recent large-scale MC simulations using 
special-purpose algorithms such as the WL technique have allowed us to settle several 
long-standing controversial questions [291 130| [33j 134"] . 
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